Genome-wide association mapping and genomic prediction for pre‑harvest sprouting resistance, low α-amylase and seed color in Iranian bread wheat

Background Pre-harvest sprouting (PHS) refers to a phenomenon, in which the physiologically mature seeds are germinated on the spike before or during the harvesting practice owing to high humidity or prolonged period of rainfall. Pre-harvest sprouting (PHS) remarkably decreases seed quality and yield in wheat; hence it is imperative to uncover genomic regions responsible for PHS tolerance to be used in wheat breeding. A genome-wide association study (GWAS) was carried out using 298 bread wheat landraces and varieties from Iran to dissect the genomic regions of PHS tolerance in a well-irrigated environment. Three different approaches (RRBLUP, GBLUP and BRR) were followed to estimate prediction accuracies in wheat genomic selection. Results Genomes B, A, and D harbored the largest number of significant marker pairs (MPs) in both landraces (427,017, 328,006, 92,702 MPs) and varieties (370,359, 266,708, 63,924 MPs), respectively. However, the LD levels were found the opposite, i.e., genomes D, A, and B have the highest LD, respectively. Association mapping by using GLM and MLM models resulted in 572 and 598 marker-trait associations (MTAs) for imputed SNPs (− log10 P > 3), respectively. Gene ontology exhibited that the pleitropic MPs located on 1A control seed color, α-Amy activity, and PHS. RRBLUP model indicated genetic effects better than GBLUP and BRR, offering a favorable tool for wheat genomic selection. Conclusions Gene ontology exhibited that the pleitropic MPs located on 1A can control seed color, α-Amy activity, and PHS. The verified markers in the current work can provide an opportunity to clone the underlying QTLs/genes, fine mapping, and genome-assisted selection.Our observations uncovered key MTAs related to seed color, α-Amy activity, and PHS that can be exploited in the genome-mediated development of novel varieties in wheat. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03628-3.

sprouting (PHS) is known as a detrimental restricting factor in wheat productivity [4]. Given this challenge, genetic improvements in PHS tolerance have become a serious focus of wheat breeders.
PHS tolerance depends on several factors, including i) environmental factors, such as relative humidity and temperature [4]; ii) biophysiological traits, such as germination-inhibitory compounds in the glumes, α-amylase (α-Amy) activity, grain structure and color, phytohormones, and seed dormancy [5]; iii) morphological traits, such as awn and spike structure [6]. Of these factors, grain color is genetically related to PHS tolerance, the red-grained genotypes are more tolerant to PHS than white ones [7]. Genes coding MYB transcriptional factors responsible for the flavonoid biosynthesis, i.e., Tamyb10-1, have been reported as candidates that determine grain color [8]. Myb10 confers PHS resistance in wheat, which activates 9-cis-epoxycarotenoid dioxygenase (NCED) by biding the secondary wall MYBresponsive element (SMRE) to promote ABA biosynthesis in early wheat seed development stages [9][10][11]. Moreover, experimental evidence highlight seed dormancy is a key genetic component that determines PHS tolerance in wheat genotypes [2].
Genome-wide association study (GWAS) is an alternative tool to determine QTLs in natural populations [15]. The establishment of genotyping technologies, from SSRs to SNPs, could facilitate association studies for accurate and efficient exploring of potential loci involved in complex traits, including PHS resistance in wheat [7,13,31] and grain-associated traits [32,33]. However, the molecular mechanisms of PHS resistance remain unclear. Genomic selection (GS) along with GWAS can dramatically accelerate genetic gain in breeding [34,35]. Several methods, including SNP-BLUP, have been suggested for genomic prediction [36].
In this study, a total of 298 Iranian wheat genotypes were evaluated for genotyping-by-sequencing (GBS)based GWAS to achieve two objectives: i) uncovering genetic loci associated with PHS resistance; (2) identifying the best model for estimating prediction accuracies in genomic selection.

Phenotypic data summary
The results of descriptive statistics of traits related to pre-harvest sprouting are shown in Table 1. Germination percentage occurred among Iranian wheat cultivars and landraces were ranged from zero to %100. The averages of germination percentage in landraces and cultivars were 71.31% and are 79.67%, respectively, which shows that native populations harbor more value of this trait. Sprouting index, sprouting score, and sprouting spike also confirm the lower pre-harvest sprouting rate of native populations than cultivated varieties. The α-Amy enzyme activities in native populations and cultivars were 9.38 and 10.76, respectively, which indicates less activity of the enzyme in landraces than that of varieties. Color indices including L, a, and b do not differ significantly between cultivars and landraces.

Assessment of SNPs
After eight Ion Proton runs, a total of 566,439,207 reads were identified with 458,363,607 (about 81%) highquality barcoded reads. A total of 133,039 unique SNPs were called after filtering out duplicated reads. After imputation and discarding the SNPs with > 20% missing values, > 10% heterozygosity, and < 5% miner allele frequency, 43,525 SNPs were identified across all 21 wheat chromosomes. Out of them, 15,951, 21,864, and 5,710 SNPs were mapped to A, B, and D genomes, respectively, which included 36.7%, 50.2%, and 13.1% of total SNPs (Fig. 2). The highest and lowest numbers of SNPs were located on 3A (4034 SNPs) and 4D (270 SNPs), respectively.

Population structure and kinship matrix
In order to determine the appropriate number of subpopulations, the number of clusters was plotted (K) against ΔK. The largest ΔK value was observed at K = 3 suggesting the presence of three subpopulations (Fig. 3a). Using the structure software, the population of 298 accessions was structured into three subpopulations, Sub1, Sub2, and Sub_3 (Fig. 2). Sub_1 contains 113 accessions with 107 landraces and 6 varieties, Sub_2 contains 111 accessions with 97 landraces and 14 varieties; Sub_3 contains 74 studies with 70 varieties and 4 landraces (Fig. 3b). Molecular markers-based PCA showed that the first and second components justified 16.9% and 6.3% of total genetic variance occurred between wheat accessions. Thus, our study can distinguish favorably cultivars and native populations (Fig. 4). As expected, a population structure was identified in the Iranian wheat landraces, with the first five eigenvalues accounting for 30.5% of genetic diversity. From the clustering results, the native populations were divided into two subgroups. Clustering based on the nearest neighbor also indicated that cultivars and landraces were appropriately separated by using the imputed markers (Fig. 5).

Linkage disequilibrium (LD)
The levels of LD in genomes A, B, and D were 2279, 1707, and 5135, respectively. This reflects that genomes D, A, and B have the highest LD, respectively (Fig. 6). An analysis on landraces identified a total of 1,867,575 marker pairs with r 2 = 0.182, of which 847,725 (45.39%) harbored significant linkages at P < 0.001. Similar to cultivars, marker pairs on chromosome 4A showed the strongest LD (r 2 = 0.369). Moreover, most of the significant marker pairs were found at distance of < 10 cM. Genomes D and B possessed the lowest and highest number of marker pairs (92,702 and 427,017), respectively. A total of 1,858,425 marker pairs with r 2 = 0.211 were identified in cultivars, of which 700,991 (37.72%) harbored significant linkages at P < 0.001. Based on the observations, most of the significant marker pairs were found at distance of < 10 cM. Genomes D and B possessed the lowest and

MTAs for morphometric seed traits
In total, 566 and 598 significant marker pairs (MTAs) were identified by using GLM and MLM approaches, respectively, for PHS-related traits (-log10 P > 3). Of the total number of MTAs in the GLM method, 204,  b, Hue, Chroma, and WI traits using the GLM method were 60,65,72,120,40,30,50,35,39,35, and 20, as well as using the MLM method were 65,66,64,170,34,30,41,35,36,37, and 20, respectively. The highest and lowest numbers of significant marker pairs using GLM and MLM methods were related to SSp (120 and 170 marker pairs) and WI (20 and 20 marker pairs), respectively. The most significant markers for PHS were on genome B, which has a greater effect on seed dormancy when compared to other genomes. However, the seed brightness (L and WI)-associated markers were located on genome A (Fig. 7). Manhattan diagrams for common areas associated with each seed trait are shown in Fig. 8.

Gene ontology
The markers with the highest significance (P < 0.0001) and pleiotropy were studied in more detail. A total of 41 markers with high significance and pleitropic were identified, most of which were on 1A, 1B, 2A, 3B, 6D, and 7A. The marker pairs located on 1A were found to be able to control seed color, α-Amy activity, and germination percentage. Some of the significant MTAs were responsible for important molecular and biological processes, including protein kinase, G protein-coupled receptor signaling, signal transduction, intracellular transport, oxidoreductase activity, Fe ion binding, oxidation-reduction process, monooxygenase activity, protein binding, regulation of transcription, and double-stranded DNA binding (Table 3).
Based on the rice reference genome, the following pathways were discovered: hormone signal transduction (

Genomic prediction
BRR, RR-BLUP, and GBLUP models using imputed SNPs exhibited the highest prediction accuracy for phenotypes 6, 3, and 2. The highest prediction accuracy by the GBLUP was achieved for SSp, Hue, and WI; by the RR-BLUP method for SS, SI, A.amy, a, L, and b; as well as by the BRR for GP and L traits (Fig. 10). BRR, RR-BLUP, and GBLUP models using significant SNPs indicated the highest prediction accuracy for phenotypes 2, 7, and 2. The highest prediction accuracy by the GBLUP was achieved for L and WI; by the RR-BLUP method for GP, SS, SI, SSp, Hue, a, and b; as well as by the BRR for A.amy trait. Overall, the RR-BLUP showed higher prediction accuracy and the BRR had a slight difference in accuracy with the RR-BLUP.

Discussion
PHS tolerance in wheat is a complicated quantitative trait influenced by genetic background and environment [4]. Thus, reliable phenotyping and genotyping for monitoring PHS tolerance can enhance the accuracy of QTL mapping. In this study, a total of 298 Iranian wheat accessions including 208 landraces and 90 cultivars were assembled as a natural population for mapping QTLs related to α-Amy enzyme activity, seed color, PHS using GWAS.
A high level of variation was uncovered in the studied traits for Iranian wheat accessions, suggesting the potential of the GWAS technique for exploring QTLs, as reported by Rahimi et al. [40]. The α-Amy enzyme activity was lower in native populations than that of cultivars. Moreover, the seeds of landraces were exposed to longer dormancy when compared to cultivars. From correlation analysis, the below facts were concluded based on the associations among α-Amy activity, grain color, and pr-harvest sprouting, i) the lower the α-Amy activity, the higher the resistance of accessions to PHS; ii) the darker the seed, the more dormant it is; and iii) the more dormant the seed, the more PHS resistant it is. Similar results were reported by Zhou et al. [3], Zhou et al. [4], and Albrecht et al. [16].
The possibility for false associations can be raised in mapping studies if population structure is not suitably accounted for [41]. Two kinds of kinships lead to a high rate of false positives in GWAS: cryptic relatedness and ancestry difference. Cryptic relatedness appears when some plant accessions are closely related; however, these shared ancestries are undisclosed to breeders [42]. Large populations inevitably consist of accessions having common ancestry from various populations. Ancestry difference also refers to various ancestries among accessions in research [33]. To evaluate the population structure in Iranian wheat accessions, PCA analysis and clustering were performed. Of results, the panel of accessions was stratified into three groups. The selection effects in breeding programs are considered as the reasons for such a genetic separation [43]. Rahimi et al. [40] observed the same grouping on these Iranian wheat accessions. Cultivars made up one group, while landraces made up the other two groups, regardless of their geographic origins. This mixture can be originated from grain exchanges between farmers in different local markets throughout the country [44].
Of the results, the detected SNPs could cover the wheat genome well. The SNPs were higher in genome B and lower in genome D. Therefore, it seems there is a direct correlation between chromosome size and SNP density [45], because of the smaller size of chromosomes B compared to A ones. The higher frequency of SNPs in genome B resulted from the evolutionary processes. This inference was also stated by Alipour et al. [46] and Mourad et al. [47].  Genomes D, A, and B have the highest LD, respectively. The strongest LD was recorded between marker pairs on chromosome 4A. The fact that cultivars exhibited higher LD in contrast to landraces, particularly in genome D, is presumably a consequence of selection throughout the time of breeding efforts for PHS-related traits [16]. The differences in LD occurred between genomes and accessions, in addition to the evolutionary processes, indicate the impact of breeding schedules. Similarly, Liu et al. [48] observed that the distance of LD decay in the native populations is less than cultivated varieties in wheat Pakistan/China collections.
Of the results, 1A, 2A, 4A, 1B, 2B, 6B, 4B, 3B, 5B, 7B, 6D, 5D, 4D, and 2D harbor genomic regions controlling PHS-related traits. Genome B possessed the highest number of MTAs, suggesting the potential of this genome in wheat adaptability to PHS. The most significant markers for PHS were on genome B, which has a greater effect on seed dormancy when compared to other genomes. However, the seed brightness-associated markers were located on genome A. These observations are in agreement with previous studies. For instance, Zhu et al. [3] mapped three key loci for PHS tolerance on chromosomes 6BL, 3BS, and 1AL, as well as validated one dCAPS and two CAPS markers for implementation in wheat genomics-based selection.
Genomic regions controlling PHS were detected in most wheat chromosomes in this study. To date, seven     [3]. Our observations showed that the darker the seed, the more dormant it will be and thus the more resistant it will be to PHS. Of justifying the cause, some associations were observed between grain color and PHS tolerance. Zhu et al. [3], for instance, discovered the positive correlations between PHS tolerance and seed color and suggested that this association occurs because the red-colored populations harbor more tolerant Qphs.ahau-1A and Qphs. ahau-3B alleles. Therefore, the authors stated that wheat seed color may be modulated collectively via Tamyb10-1 and other QTLs. In this work, MTAs related to grain color were found on 7B, 2A, etc. In this regard, the Psy1 gene coding phytoene synthase 1, responsible for yellow pigment, is co-segregated with seed brightness on 7B [49]. A major QTL controlling both a* (redness) and L* (brightness) was also reported on 2A [44]. Therefore, it seems that QTLs located on 7B and 2A are involved in wheat seed brightness, and thereby PHS tolerance. In this work, MTAs related to seed dormancy were found in some chromosomes, such as 4A. Similarly, Torada et al. [21] mapped TaMKK3-A as a candidate gene for the wheat seed dormancy, namely Phs1, on chromosomes 4A. They suggested that a single amino acid substitution in the kinase domain of this protein is related to the length of seed dormancy. From our findings, α-Amyrelated genomic regions were found on 6B, 6D, 7B, etc. This is in line with previous studies. Lazarus et al. [50] demonstrated that α-Amy-related genomic regions are multigene families located on the chromosomes 7A, 7B, 7D (α-Amy2) and 6A, 6B, 6D (α-Amy1).  TTT GGC AAT AGA TTT  GTT ATA AAC TTG ACA  ATG GCT AAG AAG  CCT CCGT   GP, SS, SI and SSp  2A  59,228  ---41 rs63525 TGC AGT TCG TAA  GCA GAG CGG CAA  TAT ACG ATA TAC CAC  TAG TAT ACT GTG TCA  CCA CTG GGGT GP, SS, SI and A.amy 7B 72,800 --- The flanking sequences of imputed SNPs were searched and aligned versus the RefSeq v1.0 ( [51], https:// urgi. versa illes. inra. fr/ blast_ iwgsc/). Interestingly, output indicated that most marker pairs are in the protein-coding regions, which control the transcription process. DNAbinding, transcription factor activity, and transmembrane transport are other examples that are likely responsible for PHS tolerance. These findings are similar to the earlier researches [31]. Based on the rice reference genome, the following pathways were discovered: metabolic pathways, hormone signal transduction, MAPK signaling pathway, purine metabolism, spliceosome, and glycolysis/gluconeogenesis. Liu et al. [52] observed that the slowed glycolysis leads to down-regulate glycerate-3-phosphate and inhibits seed germination (i.e., PHS). Torada et al. [21] uncovered a MKK3 by a map-based approach as a candidate gene for the locus Phs1 on 4A in wheat. Liu et al. [53] revealed that water status changes transcript levels of key genes involved in auxin, JA, and ethylene biosynthesis and their metabolic pathways, suggesting roles in regulating seed dormancy and germination. Nonogaki et al. [54] showed that seed germination and dormancy, the two main factors around PHS, are controled by endogenous hormone balance, especially between GA and ABA, reflecting their vital roles in PHS. Wang et al. [38] indicated that MAPK signaling and hormone signal transduction are associated with PHS. Zhang et al. [55] also highlighted that transcripts of spliceosome-related genes are abundant in the early stage of seed germination, suggesting the role of spliceosome in PHS process.
The highest prediction accuracy by GBLUP was achieved for SSp, Hue, and WI; by RR-BLUP method for SS, SI, A.amy, a, L, and b; as well as by BRR for GP and L traits. Shabannejad et al. [56] revealed BRR and RR-BLUP are superior to other GP models, which are utilized in well-irrigated and rain-fed environments, respectively. Overall, obtaining the highest GP accuracy is depend on the genomic selection method, level of LD, genetic architecture, and genetic variation [57]. In this study, RRBLUP model indicated genetic effects better than GBLUP and BRR, offering a favorable tool for wheat genomic selection. It was reported that high genetic variation would be achieved by the GBLUP if markers were closely associated with the trait of interest and/or plant populations were advanced. RR-BLUP can work well for genetic architecture consisting of numerous loci with small impacts. BRR is similar to RR-BLUP however its shrinkage depends on the size of the studied population [58].

Conclusion
In the current study, GWAS for PHS in Iranian bread wheat accessions revealed the lowest LD decay distance and the highest number of marker pairs on genome B due to evolutionary events. The loci controlling the traits of interest were mapped on 1A, 2A, 4A, 1B, 2B, 7B, 3B, 5B, 6B, 4B, 6D, 2D, 5D, and 4D. Gene ontology exhibited that the pleitropic MPs located on 1A can control seed color, α-Amy activity, and PHS. The verified markers in the current work can provide an opportunity to clone the underlying QTLs/genes, fine mapping, and genome-assisted selection.

Plant material and field trial
To monitor PHS resistance, 208 native landraces and 90 cultivars were cultured in an alpha-lattice with two repeats in three crop seasons (2017-18, 2018-19, and 2019-20) under well-irrigated conditions (
Where n represents the number of clusters, Where ni and N are the numbers of germinated and total seeds, respectively, Where mi is the number of sprouted spikes and M is the total number of spikes.
To estimate α-Amy activity, the spikes of all accessions were taken out of the refrigerator, threshing was conducted by hand to avoid damaging the seed coat or embryos. Therefore, seeds were imbibed in a petri dish for a duration of 24 h at 25 °C and then prepared for enzyme extraction [61]. 0.5 ml of the seed extract (60 mM phosphate buffer (pH 8.6) and 0.5 ml of starch solution were incubated at 37 °C for 30 min. The reaction was ceased by adding 1 ml of hydrochloric acid (0.1 N), and then 1 ml of the iodine reagent was added to the solution. The color absorption was recorded using a plate reader at 620 nm [16].

Evaluation of seed color with digital images
The digital images of wheat grains in the current work were provided by a camera (Canon SX540 HS) equipped with 800 dpi resolution. The captured images were analyzed and processed via Python 3.7 [62]. For calibration, the regression between L, a, and b indices calculated with the Japanese CR_400 colorimeter and a photo box of 17 standard colors printed on 8 cm squares were used. Where, L, a, and b are color indices.

GBS and imputation
The GBS libraries were established and sequenced for the Iranian wheat genotypes following the procedure as explained by Alipour et al. [46]. SNPs were discovered via internal alignments after trimming reads to 64 bp and categorizing them into tags. SNP calling was carried out using the UNEAK GBS pipeline, where SNPs with low allele frequency < 1% and reads with a low-quality score < 15 were removed to keep away from false-positive outputs. The imputation was accomplished according to available allele frequencies in BEAGLE version 3.3.2 [63]. The distance of LD decay was determined by the ggplot2 package in RStudio [64]. The W7984 reference genome was used because it fulfills the highest accuracy of imputation among various wheat reference genomes [65].

Population structure and kinship matrix
Population structure in the Iranian wheat accessions was assayed via STRU CTU RE version 2.3.4. An admixture model was exploited along with a simulation phase consisting of 10,000 steps for K = 1-10. In this study, ΔK was exerted to estimate the most likely number of subpopulations [66]. To measure LD among markers, the expected and observed allele frequencies were introduced into TASSEL. To determine the relationships among the Iranian wheat accessions, a neighbor-joining tree was constructed according to a pairwise distance matrix by TASSEL version 5 [67].

Genome-wide association study
The general linear model (GLM) and mixed linear model (MLM) approaches were accomplished to obtain the marker effect estimations. The GLM was performed with population structure (Q matrix) integrated as covariate to correct for the effects of population substructure. The MLM was employed with accounting for both population structure and family structure matrix (Kinship) to control both Type I and Type II errors. The association mapping was carried out using GLM and MLM functions of TASSEL5 [65,68]. To correct for multiple testing, a false discovery rate (FDR) method described was used to declare significant marker-trait associations with relevant grain phenotype descriptor. A Manhattan plot was obtained using the CMplot package to explore associations between genotypes and phenotypes.

Annotation of genes
Sequences harboring associated SNP markers were exploited for the gene annotation by aligning to the IWGSC-RefSeq V1.0 (IWGSC) using Gramene (http:// www. grame ne. org/), an integrated database for comparative genomics in plant species. The overlapping genes with the highest blast score were picked out for further analysis. The ensemble-gramene database was utilized to extract the molecular functions and biological processes of genes in the gene ontology. Moreover, the significant SNPs were utilized in the enrichment analysis of gene ontology via KOBAS version 2.0 for testing in the KEGG (https:// www. genome. jp/ kegg/).

Genomic prediction strategies
GP was calculated by various approaches, including BRR [69], GBLUP [70], and RR-BLUP [71] based on whole 43,525 marker set and GWAS on the training set. All of the analyses were performed by iPat Tool [72]. The GP accuracy was determined as Pearson's correlations (r) between GEBVs and BLUPs over the validation and training sets [73].

Statistical analysis
The descriptive statistics and correlation analysis was conducted by R 4.1 using the ggplot2, dplyr, ggpubr and psych packages to reveal the distribution of wheat traits. To classify wheat genotypes, heatmap analysis was carried out in RStudio.